Table of Contents

  • 1  Cartopy
  • 2  Cartopy Installation
  • 3  Import Modules
  • 4  Cartopy Features
  • 5  Customize Plot
  • 6  Projections
    • 6.1  Mercator
    • 6.2  PlateCarree
    • 6.3  Orthographic
    • 6.4  Robinson
    • 6.5  Help on ccrs
      • 6.5.1  Ocean modeling
    • 6.6  Review some common projections
  • 7  Regional map
    • 7.1  AlbersEqualArea projection over ethiopia
    • 7.2  Mercator projection over ethiopia
    • 7.3  PlateCarree projection over ethiopia
    • 7.4  GSHHSFeature
    • 7.5  Regional maps
    • 7.6  Plot East Africa Region
    • 7.7  Create a regional map centered over Addis Ababa
    • 7.8  Plotting a single point
    • 7.9  Plot Multiple Points
    • 7.10  Lable Data Points
  • 8  Map Tiles
    • 8.1  Zoom to specific region
  • 9  Plot netCDF
    • 9.1  Add Colorbar
    • 9.2  Add Gridlines
  • 10  Xarray with Cartopy
    • 10.1  Sub plot_1
    • 10.2  Sub plot 2
    • 10.3  Color option

Cartopy¶

  • Is a library providing cartographic tools for python.

  • It is used making publication quality maps.

  • Cartopy makes use of the powerful packages:

  • PROJ.4

  • numpy

  • shapely libraries

  • Matplotlib

  • Key features of cartopy are

    • its object oriented projection definitions.

    • its ability to transform points, lines, vectors, polygons and images between those projections.

    • over 30 different map projections.

  • Support for Basemap is expected to wrap up in 2020

Cartopy Installation¶

pip install cartopy

conda install -c conda-forge cartopy

Import Modules¶

In [1]:
import matplotlib.pylab as plt
from matplotlib import cm
import numpy as np
In [2]:
import cartopy.crs as ccrs
import cartopy.feature as cfeature

The generation of maps using cartopy can be done:

  • Create Matplotlib figure using plt.figure().

  • Add subplot to figure with projection attribute

  • Add features to map like Land, Ocean, Coastline, Borders, etc using cartopy.feature module.

  • Add your functionalities like points, lines, & surfacesetc using common matplotlib API.

  • Show plot using plt.show().

In [3]:
fig = plt.figure(figsize=(8,6))

ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())

ax.add_feature(cfeature.COASTLINE)

plt.show()

Cartopy Features¶

Adding

  • Coastline,
  • Land,
  • Lakes,
  • Borders,
  • Ocean,
  • Rivers
  • We add them to our subplot via the add_feature method.

  • We can set parameters lise linewidth, edgecolor, alpha, color, etc

| Function | Definition | |-------------------- |--------------------------------------------------- | | cfeature.BORDERS | Country boundaries | | cfeature.COASTLINE | Coastline, including major islands | | cfeature.LAKES | Natural and artificial lakes | | cfeature.LAND | Land polygons, including major islands | | cfeature.OCEAN | Ocean polygons | | cfeature.RIVERS | Single-line drainages, including lake centerlines | | cfeature.STATES | (limited to the United States at this scale) |

In [4]:
fig = plt.figure(figsize=(8,8))

ax = fig.add_subplot(1,1,1, projection=ccrs.Mercator())

ax.add_feature(cfeature.COASTLINE)

ax.add_feature(cfeature.LAND, color="lightgrey", alpha=0.5)

ax.add_feature(cfeature.LAKES, color="lime")

ax.add_feature(cfeature.BORDERS, linestyle="--")

ax.add_feature(cfeature.OCEAN, color="skyblue", alpha=0.4)

ax.add_feature(cfeature.RIVERS, edgecolor="red")

ax.add_feature(cfeature.STATES)

plt.show()

Customize Plot¶

Plot can be adjusted by:

  • set_global: zoom the map out as much as possible

  • set_extent: zoom the map to the given bounding box

  • gridlines : add a graticule (and optionally labels) to the axes

  • coastlines: add Natural Earth coastlines to the axes

  • stock_img : add a low-resolution Natural Earth background image to the axes

  • imshow : add an image (numpy array) to the axes

  • add_geometries: add a collection of geometries (Shapely) to the axes

Projections¶

Review some common projections:

  • PlateCarree

  • Robinson

  • Mercator

  • Orthographic

  • InterruptedGoodeHomolosine

Mercator¶

In [5]:
plt.figure(figsize=(6,6))
ax = plt.axes(projection=ccrs.Mercator())
ax.coastlines()
Out[5]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f93352dc730>

PlateCarree¶

In [6]:
plt.figure(figsize=(6,6))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
Out[6]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f9331554d30>

Orthographic¶

In [7]:
plt.figure(figsize=(6,6))
ax = plt.axes(projection=ccrs.Orthographic())
ax.coastlines()
Out[7]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f933157a290>

Robinson¶

In [8]:
plt.figure(figsize=(6,6))
ax = plt.axes(projection=ccrs.Robinson())
ax.coastlines()
Out[8]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f9331718400>

Help on ccrs¶

In [9]:
ccrs.PlateCarree?
In [10]:
help(ccrs.PlateCarree()) 
Help on PlateCarree in module cartopy.crs object:

class PlateCarree(_CylindricalProjection)
 |  PlateCarree(central_longitude=0.0, globe=None)
 |  
 |  Method resolution order:
 |      PlateCarree
 |      _CylindricalProjection
 |      _RectangularProjection
 |      Projection
 |      CRS
 |      pyproj.crs.crs.CustomConstructorCRS
 |      pyproj.crs.crs.CRS
 |      builtins.object
 |  
 |  Methods defined here:
 |  
 |  __init__(self, central_longitude=0.0, globe=None)
 |      Parameters
 |      ----------
 |      proj4_params: iterable of key-value pairs
 |          The proj4 parameters required to define the
 |          desired CRS.  The parameters should not describe
 |          the desired elliptic model, instead create an
 |          appropriate Globe instance. The ``proj4_params``
 |          parameters will override any parameters that the
 |          Globe defines.
 |      globe: :class:`~cartopy.crs.Globe` instance, optional
 |          If omitted, the default Globe instance will be created.
 |          See :class:`~cartopy.crs.Globe` for details.
 |  
 |  quick_vertices_transform(self, vertices, src_crs)
 |      Where possible, return a vertices array transformed to this CRS from
 |      the given vertices array of shape ``(n, 2)`` and the source CRS.
 |      
 |      Note
 |      ----
 |          This method may return None to indicate that the vertices cannot
 |          be transformed quickly, and a more complex geometry transformation
 |          is required (see :meth:`cartopy.crs.Projection.project_geometry`).
 |  
 |  ----------------------------------------------------------------------
 |  Data and other attributes defined here:
 |  
 |  __abstractmethods__ = frozenset()
 |  
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from _RectangularProjection:
 |  
 |  boundary
 |  
 |  x_limits
 |  
 |  y_limits
 |  
 |  ----------------------------------------------------------------------
 |  Methods inherited from Projection:
 |  
 |  is_geodetic(self)
 |  
 |  project_geometry(self, geometry, src_crs=None)
 |      Project the given geometry into this projection.
 |      
 |      Parameters
 |      ----------
 |      geometry
 |          The geometry to (re-)project.
 |      src_crs: optional
 |          The source CRS.  Defaults to None.
 |      
 |          If src_crs is None, the source CRS is assumed to be a geodetic
 |          version of the target CRS.
 |      
 |      Returns
 |      -------
 |      geometry
 |          The projected result (a shapely geometry).
 |  
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from Projection:
 |  
 |  ccw_boundary
 |  
 |  cw_boundary
 |  
 |  domain
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from Projection:
 |  
 |  threshold
 |  
 |  ----------------------------------------------------------------------
 |  Methods inherited from CRS:
 |  
 |  __eq__(self, other)
 |      Return self==value.
 |  
 |  __getstate__(self)
 |      Return the full state of this instance for reconstruction
 |      in ``__setstate__``.
 |  
 |  __hash__(self)
 |      Hash the CRS based on its proj4_init string.
 |  
 |  __reduce__(self)
 |      Implement the __reduce__ API so that unpickling produces a stateless
 |      instance of this class (e.g. an empty tuple). The state will then be
 |      added via __getstate__ and __setstate__.
 |      We are forced to this approach because a CRS does not store
 |      the constructor keyword arguments in its state.
 |  
 |  __setstate__(self, state)
 |      Take the dictionary created by ``__getstate__`` and passes it
 |      through to this implementation's __init__ method.
 |  
 |  as_geocentric(self)
 |      Return a new Geocentric CRS with the same ellipse/datum as this
 |      CRS.
 |  
 |  as_geodetic(self)
 |      Return a new Geodetic CRS with the same ellipse/datum as this
 |      CRS.
 |  
 |  transform_point(self, x, y, src_crs, trap=True)
 |      transform_point(x, y, src_crs)
 |      
 |      Transform the given float64 coordinate pair, in the given source
 |      coordinate system (``src_crs``), to this coordinate system.
 |      
 |      Parameters
 |      ----------
 |      x
 |          the x coordinate, in ``src_crs`` coordinates, to transform
 |      y
 |          the y coordinate, in ``src_crs`` coordinates, to transform
 |      src_crs
 |          instance of :class:`CRS` that represents the coordinate
 |          system of ``x`` and ``y``.
 |      trap
 |          Whether proj errors for "latitude or longitude exceeded limits" and
 |          "tolerance condition error" should be trapped.
 |      
 |      Returns
 |      -------
 |      (x, y) in this coordinate system
 |  
 |  transform_points(self, src_crs, x, y, z=None, trap=False)
 |      transform_points(src_crs, x, y[, z])
 |      
 |      Transform the given coordinates, in the given source
 |      coordinate system (``src_crs``), to this coordinate system.
 |      
 |      Parameters
 |      ----------
 |      src_crs
 |          instance of :class:`CRS` that represents the
 |          coordinate system of ``x``, ``y`` and ``z``.
 |      x
 |          the x coordinates (array), in ``src_crs`` coordinates,
 |          to transform.  May be 1 or 2 dimensional.
 |      y
 |          the y coordinates (array), in ``src_crs`` coordinates,
 |          to transform.  Its shape must match that of x.
 |      z: optional
 |          the z coordinates (array), in ``src_crs`` coordinates, to
 |          transform.  Defaults to None.
 |          If supplied, its shape must match that of x.
 |      trap
 |          Whether proj errors for "latitude or longitude exceeded limits" and
 |          "tolerance condition error" should be trapped.
 |      
 |      Returns
 |      -------
 |          Array of shape ``x.shape + (3, )`` in this coordinate system.
 |  
 |  transform_vectors(self, src_proj, x, y, u, v)
 |      transform_vectors(src_proj, x, y, u, v)
 |      
 |      Transform the given vector components, with coordinates in the
 |      given source coordinate system (``src_proj``), to this coordinate
 |      system. The vector components must be given relative to the
 |      source projection's coordinate reference system (grid eastward and
 |      grid northward).
 |      
 |      Parameters
 |      ----------
 |      src_proj
 |          The :class:`CRS.Projection` that represents the coordinate system
 |          the vectors are defined in.
 |      x
 |          The x coordinates of the vectors in the source projection.
 |      y
 |          The y coordinates of the vectors in the source projection.
 |      u
 |          The grid-eastward components of the vectors.
 |      v
 |          The grid-northward components of the vectors.
 |      
 |      Note
 |      ----
 |          x, y, u and v may be 1 or 2 dimensional, but must all have matching
 |          shapes.
 |      
 |      Returns
 |      -------
 |          ut, vt: The transformed vector components.
 |      
 |      Note
 |      ----
 |         The algorithm used to transform vectors is an approximation
 |         rather than an exact transform, but the accuracy should be
 |         good enough for visualization purposes.
 |  
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from CRS:
 |  
 |  proj4_params
 |  
 |  ----------------------------------------------------------------------
 |  Methods inherited from pyproj.crs.crs.CustomConstructorCRS:
 |  
 |  to_3d(self, name: Optional[str] = None) -> 'CRS'
 |      .. versionadded:: 3.1.0
 |      
 |      Convert the current CRS to the 3D version if it makes sense.
 |      
 |      New vertical axis attributes:
 |        - ellipsoidal height
 |        - oriented upwards
 |        - metre units
 |      
 |      Parameters
 |      ----------
 |      name: str, optional
 |          CRS name. Defaults to use the name of the original CRS.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  ----------------------------------------------------------------------
 |  Class methods inherited from pyproj.crs.crs.CustomConstructorCRS:
 |  
 |  from_user_input(value: Any, **kwargs) -> 'CRS' from abc.ABCMeta
 |      Initialize a CRS class instance with:
 |        - PROJ string
 |        - Dictionary of PROJ parameters
 |        - PROJ keyword arguments for parameters
 |        - JSON string with PROJ parameters
 |        - CRS WKT string
 |        - An authority string [i.e. 'epsg:4326']
 |        - An EPSG integer code [i.e. 4326]
 |        - A tuple of ("auth_name": "auth_code") [i.e ('epsg', '4326')]
 |        - An object with a `to_wkt` method.
 |        - A :class:`pyproj.crs.CRS` class
 |      
 |      Parameters
 |      ----------
 |      value : obj
 |          A Python int, dict, or str.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from pyproj.crs.crs.CustomConstructorCRS:
 |  
 |  geodetic_crs
 |      .. versionadded:: 2.2.0
 |      
 |      Returns
 |      -------
 |      CRS:
 |          The geodeticCRS / geographicCRS from the CRS.
 |  
 |  source_crs
 |      The base CRS of a BoundCRS or a DerivedCRS/ProjectedCRS,
 |      or the source CRS of a CoordinateOperation.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  sub_crs_list
 |      If the CRS is a compound CRS, it will return a list of sub CRS objects.
 |      
 |      Returns
 |      -------
 |      List[CRS]
 |  
 |  target_crs
 |      .. versionadded:: 2.2.0
 |      
 |      Returns
 |      -------
 |      CRS:
 |          The hub CRS of a BoundCRS or the target CRS of a CoordinateOperation.
 |  
 |  ----------------------------------------------------------------------
 |  Methods inherited from pyproj.crs.crs.CRS:
 |  
 |  __repr__(self) -> str
 |      Return repr(self).
 |  
 |  __str__(self) -> str
 |      Return str(self).
 |  
 |  cs_to_cf(self) -> List[dict]
 |      .. versionadded:: 3.0.0
 |      
 |      This converts all coordinate systems (cs) in the CRS
 |      to a list of Climate and Forecast (CF) Version 1.8 dicts.
 |      
 |      :ref:`build_crs_cf`
 |      
 |      Returns
 |      -------
 |      List[dict]:
 |          CF-1.8 version of the coordinate systems.
 |  
 |  equals(self, other: Any, ignore_axis_order: bool = False) -> bool
 |      .. versionadded:: 2.5.0
 |      
 |      Check if the CRS objects are equivalent.
 |      
 |      Parameters
 |      ----------
 |      other: Any
 |          Check if the other object is equivalent to this object.
 |          If the other object is not a CRS, it will try to create one.
 |          On Failure, it will return False.
 |      ignore_axis_order: bool, default=False
 |          If True, it will compare the CRS class and ignore the axis order.
 |      
 |      Returns
 |      -------
 |      bool
 |  
 |  get_geod(self) -> Optional[pyproj.geod.Geod]
 |      Returns
 |      -------
 |      pyproj.geod.Geod:
 |          Geod object based on the ellipsoid.
 |  
 |  is_exact_same(self, other: Any) -> bool
 |      Check if the CRS objects are the exact same.
 |      
 |      Parameters
 |      ----------
 |      other: Any
 |          Check if the other CRS is the exact same to this object.
 |          If the other object is not a CRS, it will try to create one.
 |          On Failure, it will return False.
 |      
 |      Returns
 |      -------
 |      bool
 |  
 |  list_authority(self, auth_name: Optional[str] = None, min_confidence: int = 70) -> List[importlib._bootstrap.AuthorityMatchInfo]
 |      .. versionadded:: 3.2.0
 |      
 |      Return the authority names and codes best matching the CRS.
 |      
 |      Example:
 |      
 |      >>> from pyproj import CRS
 |      >>> ccs = CRS("epsg:4328")
 |      >>> ccs.list_authority()
 |      [AuthorityMatchInfo(auth_name='EPSG', code='4326', confidence=100)]
 |      
 |      If the CRS is bound, you can get an authority from
 |      the source CRS:
 |      
 |      >>> from pyproj import CRS
 |      >>> ccs = CRS("+proj=geocent +datum=WGS84 +towgs84=0,0,0")
 |      >>> ccs.list_authority()
 |      []
 |      >>> ccs.source_crs.list_authority()
 |      [AuthorityMatchInfo(auth_name='EPSG', code='4978', confidence=70)]
 |      >>> ccs == CRS.from_authorty('EPSG', '4978')
 |      False
 |      
 |      Parameters
 |      ----------
 |      auth_name: str, optional
 |          The name of the authority to filter by.
 |      min_confidence: int, default=70
 |          A value between 0-100 where 100 is the most confident.
 |          :ref:`min_confidence`
 |      
 |      Returns
 |      -------
 |      List[AuthorityMatchInfo]:
 |          List of authority matches for the CRS.
 |  
 |  to_authority(self, auth_name: Optional[str] = None, min_confidence: int = 70)
 |      .. versionadded:: 2.2.0
 |      
 |      Return the authority name and code best matching the CRS
 |      or None if it a match is not found.
 |      
 |      Example:
 |      
 |      >>> from pyproj import CRS
 |      >>> ccs = CRS("epsg:4328")
 |      >>> ccs.to_authority()
 |      ('EPSG', '4328')
 |      
 |      If the CRS is bound, you can get an authority from
 |      the source CRS:
 |      
 |      >>> from pyproj import CRS
 |      >>> ccs = CRS("+proj=geocent +datum=WGS84 +towgs84=0,0,0")
 |      >>> ccs.to_authority()
 |      >>> ccs.source_crs.to_authority()
 |      ('EPSG', '4978')
 |      >>> ccs == CRS.from_authorty('EPSG', '4978')
 |      False
 |      
 |      Parameters
 |      ----------
 |      auth_name: str, optional
 |          The name of the authority to filter by.
 |      min_confidence: int, default=70
 |          A value between 0-100 where 100 is the most confident.
 |          :ref:`min_confidence`
 |      
 |      Returns
 |      -------
 |      tuple(str, str) or None:
 |          The best matching (<auth_name>, <code>) for the confidence level.
 |  
 |  to_cf(self, wkt_version: Union[pyproj.enums.WktVersion, str] = <WktVersion.WKT2_2019: 'WKT2_2019'>, errcheck: bool = False) -> dict
 |      .. versionadded:: 2.2.0
 |      
 |      This converts a :obj:`pyproj.crs.CRS` object
 |      to a Climate and Forecast (CF) Grid Mapping Version 1.8 dict.
 |      
 |      :ref:`build_crs_cf`
 |      
 |      Parameters
 |      ----------
 |      wkt_version: str or pyproj.enums.WktVersion
 |          Version of WKT supported by CRS.to_wkt.
 |          Default is :attr:`pyproj.enums.WktVersion.WKT2_2019`.
 |      errcheck: bool, default=False
 |          If True, will warn when parameters are ignored.
 |      
 |      Returns
 |      -------
 |      dict:
 |          CF-1.8 version of the projection.
 |  
 |  to_dict(self) -> dict
 |      .. versionadded:: 2.2.0
 |      
 |      Converts the CRS to dictionary of PROJ parameters.
 |      
 |      .. warning:: You will likely lose important projection
 |        information when converting to a PROJ string from
 |        another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems  # noqa: E501
 |      
 |      Returns
 |      -------
 |      dict:
 |          PROJ params in dict format.
 |  
 |  to_epsg(self, min_confidence: int = 70) -> Optional[int]
 |      Return the EPSG code best matching the CRS
 |      or None if it a match is not found.
 |      
 |      Example:
 |      
 |      >>> from pyproj import CRS
 |      >>> ccs = CRS("epsg:4328")
 |      >>> ccs.to_epsg()
 |      4328
 |      
 |      If the CRS is bound, you can attempt to get an epsg code from
 |      the source CRS:
 |      
 |      >>> from pyproj import CRS
 |      >>> ccs = CRS("+proj=geocent +datum=WGS84 +towgs84=0,0,0")
 |      >>> ccs.to_epsg()
 |      >>> ccs.source_crs.to_epsg()
 |      4978
 |      >>> ccs == CRS.from_epsg(4978)
 |      False
 |      
 |      Parameters
 |      ----------
 |      min_confidence: int, default=70
 |          A value between 0-100 where 100 is the most confident.
 |          :ref:`min_confidence`
 |      
 |      
 |      Returns
 |      -------
 |      Optional[int]:
 |          The best matching EPSG code matching the confidence level.
 |  
 |  to_json(self, pretty: bool = False, indentation: int = 2) -> str
 |      .. versionadded:: 2.4.0
 |      
 |      Convert the object to a JSON string.
 |      
 |      Parameters
 |      ----------
 |      pretty: bool, default=False
 |          If True, it will set the output to be a multiline string.
 |      indentation: int, default=2
 |          If pretty is True, it will set the width of the indentation.
 |      
 |      Returns
 |      -------
 |      str
 |  
 |  to_json_dict(self) -> dict
 |      .. versionadded:: 2.4.0
 |      
 |      Convert the object to a JSON dictionary.
 |      
 |      Returns
 |      -------
 |      dict
 |  
 |  to_proj4(self, version: Union[pyproj.enums.ProjVersion, int] = <ProjVersion.PROJ_5: 5>) -> str
 |      Convert the projection to a PROJ string.
 |      
 |      .. warning:: You will likely lose important projection
 |        information when converting to a PROJ string from
 |        another format. See:
 |        https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems  # noqa: E501
 |      
 |      Parameters
 |      ----------
 |      version: pyproj.enums.ProjVersion
 |          The version of the PROJ string output.
 |          Default is :attr:`pyproj.enums.ProjVersion.PROJ_4`.
 |      
 |      Returns
 |      -------
 |      str
 |  
 |  to_string(self) -> str
 |      .. versionadded:: 2.2.0
 |      
 |      Convert the CRS to a string.
 |      
 |      It attempts to convert it to the authority string.
 |      Otherwise, it uses the string format of the user
 |      input to create the CRS.
 |      
 |      Returns
 |      -------
 |      str
 |  
 |  to_wkt(self, version: Union[pyproj.enums.WktVersion, str] = <WktVersion.WKT2_2019: 'WKT2_2019'>, pretty: bool = False) -> str
 |      Convert the projection to a WKT string.
 |      
 |      Version options:
 |        - WKT2_2015
 |        - WKT2_2015_SIMPLIFIED
 |        - WKT2_2019
 |        - WKT2_2019_SIMPLIFIED
 |        - WKT1_GDAL
 |        - WKT1_ESRI
 |      
 |      
 |      Parameters
 |      ----------
 |      version: pyproj.enums.WktVersion, optional
 |          The version of the WKT output.
 |          Default is :attr:`pyproj.enums.WktVersion.WKT2_2019`.
 |      pretty: bool, default=False
 |          If True, it will set the output to be a multiline string.
 |      
 |      Returns
 |      -------
 |      str
 |  
 |  ----------------------------------------------------------------------
 |  Class methods inherited from pyproj.crs.crs.CRS:
 |  
 |  from_authority(auth_name: str, code: Union[str, int]) -> 'CRS' from abc.ABCMeta
 |      .. versionadded:: 2.2.0
 |      
 |      Make a CRS from an authority name and authority code
 |      
 |      Parameters
 |      ----------
 |      auth_name: str
 |          The name of the authority.
 |      code : int or str
 |          The code used by the authority.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  from_dict(proj_dict: dict) -> 'CRS' from abc.ABCMeta
 |      .. versionadded:: 2.2.0
 |      
 |      Make a CRS from a dictionary of PROJ parameters.
 |      
 |      Parameters
 |      ----------
 |      proj_dict : str
 |          PROJ params in dict format.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  from_epsg(code: Union[str, int]) -> 'CRS' from abc.ABCMeta
 |      Make a CRS from an EPSG code
 |      
 |      Parameters
 |      ----------
 |      code : int or str
 |          An EPSG code.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  from_json(crs_json: str) -> 'CRS' from abc.ABCMeta
 |      .. versionadded:: 2.4.0
 |      
 |      Create CRS from a CRS JSON string.
 |      
 |      Parameters
 |      ----------
 |      crs_json: str
 |          CRS JSON string.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  from_json_dict(crs_dict: dict) -> 'CRS' from abc.ABCMeta
 |      .. versionadded:: 2.4.0
 |      
 |      Create CRS from a JSON dictionary.
 |      
 |      Parameters
 |      ----------
 |      crs_dict: dict
 |          CRS dictionary.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  from_proj4(in_proj_string: str) -> 'CRS' from abc.ABCMeta
 |      .. versionadded:: 2.2.0
 |      
 |      Make a CRS from a PROJ string
 |      
 |      Parameters
 |      ----------
 |      in_proj_string : str
 |          A PROJ string.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  from_string(in_crs_string: str) -> 'CRS' from abc.ABCMeta
 |      Make a CRS from:
 |      
 |      Initialize a CRS class instance with:
 |       - PROJ string
 |       - JSON string with PROJ parameters
 |       - CRS WKT string
 |       - An authority string [i.e. 'epsg:4326']
 |      
 |      Parameters
 |      ----------
 |      in_crs_string : str
 |          An EPSG, PROJ, or WKT string.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  from_wkt(in_wkt_string: str) -> 'CRS' from abc.ABCMeta
 |      .. versionadded:: 2.2.0
 |      
 |      Make a CRS from a WKT string
 |      
 |      Parameters
 |      ----------
 |      in_wkt_string : str
 |          A WKT string.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  ----------------------------------------------------------------------
 |  Static methods inherited from pyproj.crs.crs.CRS:
 |  
 |  from_cf(in_cf: dict, ellipsoidal_cs: Any = None, cartesian_cs: Any = None, vertical_cs: Any = None, errcheck=False) -> 'CRS'
 |      .. versionadded:: 2.2.0
 |      
 |      .. versionadded:: 3.0.0 ellipsoidal_cs, cartesian_cs, vertical_cs
 |      
 |      .. deprecated:: 3.2.0 errcheck
 |      
 |      This converts a Climate and Forecast (CF) Grid Mapping Version 1.8
 |      dict to a :obj:`pyproj.crs.CRS` object.
 |      
 |      :ref:`build_crs_cf`
 |      
 |      Parameters
 |      ----------
 |      in_cf: dict
 |          CF version of the projection.
 |      ellipsoidal_cs: Any, optional
 |          Input to create an Ellipsoidal Coordinate System.
 |          Anything accepted by :meth:`pyproj.crs.CoordinateSystem.from_user_input`
 |          or an Ellipsoidal Coordinate System created from :ref:`coordinate_system`.
 |      cartesian_cs: Any, optional
 |          Input to create a Cartesian Coordinate System.
 |          Anything accepted by :meth:`pyproj.crs.CoordinateSystem.from_user_input`
 |          or :class:`pyproj.crs.coordinate_system.Cartesian2DCS`.
 |      vertical_cs: Any, optional
 |          Input to create a Vertical Coordinate System accepted by
 |          :meth:`pyproj.crs.CoordinateSystem.from_user_input`
 |          or :class:`pyproj.crs.coordinate_system.VerticalCS`
 |      errcheck: bool, default=False
 |          This parameter is for backwards compatibility with the old version.
 |          It currently does nothing when True or False. DEPRECATED.
 |      
 |      Returns
 |      -------
 |      CRS
 |  
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from pyproj.crs.crs.CRS:
 |  
 |  area_of_use
 |      Returns
 |      -------
 |      AreaOfUse:
 |          The area of use object with associated attributes.
 |  
 |  axis_info
 |      Retrieves all relevant axis information in the CRS.
 |      If it is a Bound CRS, it gets the axis list from the Source CRS.
 |      If it is a Compound CRS, it gets the axis list from the Sub CRS list.
 |      
 |      Returns
 |      -------
 |      List[Axis]:
 |          The list of axis information.
 |  
 |  coordinate_operation
 |      .. versionadded:: 2.2.0
 |      
 |      Returns
 |      -------
 |      CoordinateOperation
 |  
 |  coordinate_system
 |      .. versionadded:: 2.2.0
 |      
 |      Returns
 |      -------
 |      CoordinateSystem
 |  
 |  datum
 |      .. versionadded:: 2.2.0
 |      
 |      Returns
 |      -------
 |      Datum
 |  
 |  ellipsoid
 |      .. versionadded:: 2.2.0
 |      
 |      Returns
 |      -------
 |      Ellipsoid:
 |          The ellipsoid object with associated attributes.
 |  
 |  is_bound
 |      Returns
 |      -------
 |      bool:
 |          True if CRS is bound.
 |  
 |  is_compound
 |      .. versionadded:: 3.1.0
 |      
 |      Returns
 |      -------
 |      bool:
 |          True if CRS is compound.
 |  
 |  is_derived
 |      .. versionadded:: 3.2.0
 |      
 |      Returns
 |      -------
 |      bool:
 |          True if CRS is a Derived CRS.
 |  
 |  is_engineering
 |      .. versionadded:: 2.2.0
 |      
 |      Returns
 |      -------
 |      bool:
 |          True if CRS is local/engineering.
 |  
 |  is_geocentric
 |      This checks if the CRS is geocentric and
 |      takes into account if the CRS is bound.
 |      
 |      Returns
 |      -------
 |      bool:
 |          True if CRS is in geocentric (x/y) coordinates.
 |  
 |  is_geographic
 |      This checks if the CRS is geographic.
 |      It will check if it has a geographic CRS
 |      in the sub CRS if it is a compount CRS and will check if
 |      the source CRS is geographic if it is a bound CRS.
 |      
 |      Returns
 |      -------
 |      bool:
 |          True if the CRS is in geographic (lon/lat) coordinates.
 |  
 |  is_projected
 |      This checks if the CRS is projected.
 |      It will check if it has a projected CRS
 |      in the sub CRS if it is a compount CRS and will check if
 |      the source CRS is projected if it is a bound CRS.
 |      
 |      Returns
 |      -------
 |      bool:
 |          True if CRS is projected.
 |  
 |  is_vertical
 |      .. versionadded:: 2.2.0
 |      
 |      This checks if the CRS is vertical.
 |      It will check if it has a vertical CRS
 |      in the sub CRS if it is a compount CRS and will check if
 |      the source CRS is vertical if it is a bound CRS.
 |      
 |      Returns
 |      -------
 |      bool:
 |          True if CRS is vertical.
 |  
 |  name
 |      Returns
 |      -------
 |      str:
 |          The name of the CRS (from :cpp:func:`proj_get_name`).
 |  
 |  prime_meridian
 |      .. versionadded:: 2.2.0
 |      
 |      Returns
 |      -------
 |      PrimeMeridian:
 |          The prime meridian object with associated attributes.
 |  
 |  remarks
 |      .. versionadded:: 2.4.0
 |      
 |      Returns
 |      -------
 |      str:
 |          Remarks about object.
 |  
 |  scope
 |      .. versionadded:: 2.4.0
 |      
 |      Returns
 |      -------
 |      str:
 |          Scope of object.
 |  
 |  type_name
 |      Returns
 |      -------
 |      str:
 |          The name of the type of the CRS object.
 |  
 |  utm_zone
 |      .. versionadded:: 2.6.0
 |      
 |      Finds the UTM zone in a Projected CRS, Bound CRS, or Compound CRS
 |      
 |      Returns
 |      -------
 |      Optional[str]:
 |          The UTM zone number and letter if applicable.
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from pyproj.crs.crs.CRS:
 |  
 |  __dict__
 |      dictionary for instance variables (if defined)
 |  
 |  __weakref__
 |      list of weak references to the object (if defined)

Ocean modeling¶

In [11]:
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
Out[11]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f93317184f0>
In [12]:
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
ax.coastlines()
Out[12]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7f933172b850>

Review some common projections¶

In [13]:
projections = [ccrs.PlateCarree(), ccrs.Robinson(), ccrs.Mercator(), ccrs.Orthographic(), 
               ccrs.InterruptedGoodeHomolosine()]

for proj in projections:
    plt.figure()
    ax = plt.axes(projection=proj)
    ax.stock_img()
    ax.coastlines()
    ax.gridlines()
    ax.set_global() 
    ax.set_title(f'{type(proj)}')

List of cartopy projections: https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html#

Regional map¶

  • use set_extent to limit the size of the region
In [14]:
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Mercator())
ax.coastlines()
ax.set_extent([33, 48, 3, 15])

AlbersEqualArea projection over ethiopia¶

In [15]:
plt.figure(figsize=(6, 6))

central_lon, central_lat = 37,7
extent = [33, 48, 3, 15]

ax1 = plt.axes(projection=ccrs.AlbersEqualArea(central_lon, central_lat))

ax1.set_extent(extent)

ax1.coastlines(resolution='10m')

ax1.add_feature(cfeature.OCEAN)
ax1.add_feature(cfeature.BORDERS, edgecolor='black')
ax1.add_feature(cfeature.LAND, edgecolor='black')
ax1.add_feature(cfeature.LAKES, edgecolor='black')
ax1.add_feature(cfeature.RIVERS)

ax1.gridlines()
Out[15]:
<cartopy.mpl.gridliner.Gridliner at 0x7f930fcf5720>

Mercator projection over ethiopia¶

In [16]:
plt.figure(figsize=(6, 6))

central_lon, central_lat = 37,7
extent = [33, 48, 3, 15]

ax3 = plt.axes(projection=ccrs.Mercator(central_longitude=central_lon, min_latitude=3, max_latitude=15, 
                                        globe=None, latitude_true_scale=None, false_easting=0.0, false_northing=0.0, scale_factor=None))

ax3.set_extent(extent)

ax3.coastlines(resolution='50m')

ax3.add_feature(cfeature.OCEAN)
ax3.add_feature(cfeature.BORDERS, edgecolor='black')
ax3.add_feature(cfeature.LAND, edgecolor='black')
ax3.add_feature(cfeature.LAKES, edgecolor='black')
ax3.add_feature(cfeature.RIVERS)

ax3.gridlines()
Out[16]:
<cartopy.mpl.gridliner.Gridliner at 0x7f930e4b6110>

PlateCarree projection over ethiopia¶

In [17]:
plt.figure(figsize=(6, 6))

central_lon, central_lat = 37,7
extent = [33, 48, 3, 15]

ax2 = plt.axes(projection=ccrs.PlateCarree(central_longitude=central_lon, globe=None))
ax2.set_extent(extent)

ax2.coastlines(resolution='110m')

ax2.add_feature(cfeature.OCEAN)
ax2.add_feature(cfeature.BORDERS, edgecolor='black')
ax2.add_feature(cfeature.LAND, edgecolor='black')
ax2.add_feature(cfeature.LAKES, edgecolor='black')
ax2.add_feature(cfeature.RIVERS)
ax2.gridlines()
Out[17]:
<cartopy.mpl.gridliner.Gridliner at 0x7f930e5a55a0>

GSHHSFeature¶

GSHHG (Global Self-consistent, Hierarchical, High-resolution Geography Database

In [18]:
plt.figure(figsize=(6,6))
ax = plt.axes(projection=ccrs.Mercator())
ax.coastlines(resolution = '50m')
ax.set_extent([33, 48, 3, 15], crs=ccrs.PlateCarree())

ax.add_feature(cfeature.GSHHSFeature(levels=[1], 
                                     scale="intermediate", 
                                     facecolor="lightgray"))
# scale = ‘auto’, ‘coarse’, ‘low’, ‘intermediate’, ‘high, or ‘full’ 
# levels: list of integers corresponding to the desired GSHHS feature

#ax.coastlines(resolution='50m')

ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.BORDERS, edgecolor='black')
ax.add_feature(cfeature.LAND, edgecolor='black')
ax.add_feature(cfeature.LAKES, edgecolor='black')
#ax.add_feature(cfeature.RIVERS)
ax.gridlines()
Out[18]:
<cartopy.mpl.gridliner.Gridliner at 0x7f930f5aafb0>

Regional maps¶

In [19]:
projPC = ccrs.PlateCarree()
lonW = 33
lonE = 48
latS = 3
latN = 15
cLat = (latN + latS) / 2
cLon = (lonW + lonE) / 2
res = '110m'
In [20]:
fig = plt.figure(figsize=(11, 8.5))

ax = plt.subplot(1, 1, 1, projection=projPC)

ax.set_title('Plate Carree')

gl = ax.gridlines(draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')

ax.set_extent([lonW, lonE, latS, latN], crs=projPC)
ax.coastlines(resolution=res, color='black')

ax.add_feature(cfeature.STATES, linewidth=0.3, edgecolor='brown')
ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue');

Plot East Africa Region¶

In [21]:
fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())

ax.set_extent([20, 55, -10, 20], crs=ccrs.PlateCarree())
# Put a background image on for nice sea rendering.

ax.stock_img()

ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.COASTLINE)

ax.gridlines(xlocs=np.arange(15,60,5), ylocs=np.arange(-10,20,5), draw_labels=True, crs=ccrs.PlateCarree())
Out[21]:
<cartopy.mpl.gridliner.Gridliner at 0x7f930eb6bbe0>

Create a regional map centered over Addis Ababa¶

In [22]:
latN = 9.20
latS = 8.80
lonW = 39.00
lonE = 38.55
cLat = (latN + latS) / 2
cLon = (lonW + lonE) / 2
projadd = ccrs.Mercator()
In [23]:
fig = plt.figure(figsize=(15, 10))
ax = plt.subplot(1, 1, 1, projection=projadd)
ax.set_extent([lonW, lonE, latS, latN], crs=projPC)
ax.set_facecolor(cfeature.COLORS['water'])

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle='--')
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.STATES)
ax.add_feature(cfeature.RIVERS)

gl = ax.gridlines(draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')


ax.set_title('Addis Ababa');

Plotting a single point¶

In [24]:
plt.figure(figsize=(6,6))

lat = 10
lon = 39

ax = plt.axes(projection=ccrs.Mercator())
ax.set_extent([33.,48.,3.,15.])

ax.coastlines(resolution='50m')

ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.BORDERS, edgecolor='black')
ax.add_feature(cfeature.LAND, edgecolor='black')
ax.add_feature(cfeature.LAKES, edgecolor='black')
#ax.add_feature(cfeature.RIVERS)

ax.scatter(lon,lat,50,transform=ccrs.Geodetic())
Out[24]:
<matplotlib.collections.PathCollection at 0x7f930ec38af0>

Plot Multiple Points¶

In [25]:
# Define the data points
latitudes = [10.90, 3.82, 6.60, 12.68, 8.09, 12.51, 14.21, 11.38, 4.19, 13.33]
longitudes = [34.41, 36.36, 40.55, 47.89, 44.39, 44.13, 44.9, 40.36, 39.14, 47.05]
data = [200, 150, 175, 189, 100, 315, 222, 275, 155, 190]
In [26]:
# Create a figure and axes
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# Set the extent of the map
ax.set_extent([33, 48, 3, 15], crs=ccrs.PlateCarree())

# Add map features
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle='-', alpha=0.5)
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.RIVERS)

# Create a dummy image to use for the colorbar
im = ax.imshow([[np.min(data), np.max(data)]], cmap='cool')

# Remove the dummy image from the plot
im.remove()

# Plot the data points as markers on the map
sc = ax.scatter(longitudes, latitudes, c=data, cmap='cool', s=100, transform=ccrs.PlateCarree())

# Add a colorbar
plt.colorbar(sc, label='Data')

# Add a title
plt.title('Data Points ')

# Show the plot
plt.show()

Lable Data Points¶

In [27]:
# Create a figure and axes
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# Set the extent of the map
ax.set_extent([33, 48, 3, 15], crs=ccrs.PlateCarree())

# Add map features
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle='-', alpha=0.5)
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.RIVERS)

# Create a dummy image to use for the colorbar
im = ax.imshow([[np.min(data), np.max(data)]], cmap='cool')

# Remove the dummy image from the plot
im.remove()

# Plot the data points as markers on the map
sc = ax.scatter(longitudes, latitudes, c=data, cmap='cool', s=100, transform=ccrs.PlateCarree())

# Add labels to the data points
for i in range(len(latitudes)):
    ax.text(longitudes[i], latitudes[i], str(data[i]), ha='center', va='center', transform=ccrs.PlateCarree())

# Add a colorbar
plt.colorbar(sc, label='Data')

# Add a title
plt.title('Label Data')

# Show the plot
plt.show()

Map Tiles¶

In [28]:
from cartopy.io.img_tiles import GoogleTiles, OSM,QuadtreeTiles,MapQuestOpenAerial,MapQuestOSM
In [29]:
plt.figure(figsize=(8,7))

lat = 9
lon = 39

ax = plt.axes(projection=ccrs.Mercator())
ax.set_extent([33.,48.,3.,15.])

ax.coastlines(resolution='50m',zorder=9)

ax.scatter(lon,lat,80,transform=ccrs.PlateCarree(),zorder=10)

gg_tiles = OSM()

ax.add_image(gg_tiles, 5)

Zoom to specific region¶

In [30]:
plt.figure(figsize=(6,6))

lat = 9
lon = 39

ax = plt.axes(projection=ccrs.Mercator())
ax.set_extent([38.9,39.1, 8.9,9.1])

ax.scatter(lon,lat,100,transform=ccrs.Geodetic(),zorder=10)

gg_tiles = OSM()
ax.add_image(gg_tiles, 12)

Plot netCDF¶

In [31]:
from netCDF4 import Dataset
In [32]:
fl = Dataset('era5_ecmwf_mon.nc')
In [33]:
for v in fl.variables:
    print(v)
longitude
latitude
time
t2m
sst
tp
In [34]:
print(fl.variables['t2m'])
<class 'netCDF4._netCDF4.Variable'>
int16 t2m(time, latitude, longitude)
    scale_factor: 0.0004388592865426579
    add_offset: 293.35793486723173
    _FillValue: -32767
    missing_value: -32767
    units: K
    long_name: 2 metre temperature
unlimited dimensions: 
current shape = (60, 49, 61)
filling on
In [35]:
print(fl.variables['sst'])
<class 'netCDF4._netCDF4.Variable'>
int16 sst(time, latitude, longitude)
    scale_factor: 0.0002082942408309554
    add_offset: 299.6367366731921
    _FillValue: -32767
    missing_value: -32767
    units: K
    long_name: Sea surface temperature
unlimited dimensions: 
current shape = (60, 49, 61)
filling on
In [36]:
print(fl.variables['tp'])
<class 'netCDF4._netCDF4.Variable'>
int16 tp(time, latitude, longitude)
    scale_factor: 2.9053319214141008e-08
    add_offset: 0.0009519610573705443
    _FillValue: -32767
    missing_value: -32767
    units: m
    long_name: Total precipitation
unlimited dimensions: 
current shape = (60, 49, 61)
filling on
In [37]:
print(fl.variables['time'])
<class 'netCDF4._netCDF4.Variable'>
int32 time(time)
    units: hours since 1900-01-01 00:00:00.0
    long_name: time
    calendar: gregorian
unlimited dimensions: 
current shape = (60,)
filling on, default _FillValue of -2147483647 used
In [38]:
air = fl.variables['t2m'][0,:,:]
lat = fl.variables['latitude'][:]
lon = fl.variables['longitude'][:]
In [39]:
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Mercator())
ax.coastlines(resolution = '50m')
ax.set_extent([33, 48, 3, 15], crs=ccrs.PlateCarree())

ax.coastlines(resolution='50m')

ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.BORDERS, edgecolor='black')
ax.add_feature(cfeature.LAND, edgecolor='black')
ax.add_feature(cfeature.LAKES, edgecolor='black')

ax.contourf(lon, lat, air, transform=ccrs.PlateCarree())
Out[39]:
<cartopy.mpl.contour.GeoContourSet at 0x7f930b540b50>

Add Colorbar¶

In [40]:
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Mercator())
ax.coastlines(resolution = '50m')
ax.set_extent([33, 48, 3, 15], crs=ccrs.PlateCarree())

ax.coastlines(resolution='50m')

ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.BORDERS, edgecolor='black')
ax.add_feature(cfeature.LAND, edgecolor='black')
ax.add_feature(cfeature.LAKES, edgecolor='black')

image = ax.contourf(lon, lat, air, transform=ccrs.PlateCarree())

cb = plt.colorbar(image, orientation='vertical', pad=0.03, shrink = 0.5)
cb.ax.tick_params(labelsize=12)
cb.set_label('Kelvin [K]', size=10)

Add Gridlines¶

In [41]:
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Mercator())
ax.coastlines(resolution = '50m')
ax.set_extent([33, 48, 3, 15], crs=ccrs.PlateCarree())

ax.coastlines(resolution='50m')

ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.BORDERS, edgecolor='black')
ax.add_feature(cfeature.LAND, edgecolor='black')
ax.add_feature(cfeature.LAKES, edgecolor='black')

image = ax.contourf(lon, lat, air, transform=ccrs.PlateCarree())

cb = plt.colorbar(image, orientation='horizontal', pad=0.09, shrink = 0.5)
cb.ax.tick_params(labelsize=12)
cb.set_label('Kelvin [K]', size=10)

# Only PlateCarree and Mercator plots are currently supported.
gl = ax.gridlines(draw_labels=True) 

Xarray with Cartopy¶

In [42]:
import xarray as xr
In [43]:
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter 
In [44]:
era5_ecmwf_mon_ds = xr.open_dataset('era5_ecmwf_mon.nc')
era5_ecmwf_mon_ds
Out[44]:
<xarray.Dataset>
Dimensions:    (longitude: 61, latitude: 49, time: 60)
Coordinates:
  * longitude  (longitude) float32 33.0 33.25 33.5 33.75 ... 47.5 47.75 48.0
  * latitude   (latitude) float32 15.0 14.75 14.5 14.25 ... 3.75 3.5 3.25 3.0
  * time       (time) datetime64[ns] 2018-01-01T23:00:00 ... 2022-12-01T23:00:00
Data variables:
    t2m        (time, latitude, longitude) float32 ...
    sst        (time, latitude, longitude) float32 ...
    tp         (time, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2023-03-14 04:44:29 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...
xarray.Dataset
    • longitude: 61
    • latitude: 49
    • time: 60
    • longitude
      (longitude)
      float32
      33.0 33.25 33.5 ... 47.5 47.75 48.0
      units :
      degrees_east
      long_name :
      longitude
      array([33.  , 33.25, 33.5 , 33.75, 34.  , 34.25, 34.5 , 34.75, 35.  , 35.25,
             35.5 , 35.75, 36.  , 36.25, 36.5 , 36.75, 37.  , 37.25, 37.5 , 37.75,
             38.  , 38.25, 38.5 , 38.75, 39.  , 39.25, 39.5 , 39.75, 40.  , 40.25,
             40.5 , 40.75, 41.  , 41.25, 41.5 , 41.75, 42.  , 42.25, 42.5 , 42.75,
             43.  , 43.25, 43.5 , 43.75, 44.  , 44.25, 44.5 , 44.75, 45.  , 45.25,
             45.5 , 45.75, 46.  , 46.25, 46.5 , 46.75, 47.  , 47.25, 47.5 , 47.75,
             48.  ], dtype=float32)
    • latitude
      (latitude)
      float32
      15.0 14.75 14.5 ... 3.5 3.25 3.0
      units :
      degrees_north
      long_name :
      latitude
      array([15.  , 14.75, 14.5 , 14.25, 14.  , 13.75, 13.5 , 13.25, 13.  , 12.75,
             12.5 , 12.25, 12.  , 11.75, 11.5 , 11.25, 11.  , 10.75, 10.5 , 10.25,
             10.  ,  9.75,  9.5 ,  9.25,  9.  ,  8.75,  8.5 ,  8.25,  8.  ,  7.75,
              7.5 ,  7.25,  7.  ,  6.75,  6.5 ,  6.25,  6.  ,  5.75,  5.5 ,  5.25,
              5.  ,  4.75,  4.5 ,  4.25,  4.  ,  3.75,  3.5 ,  3.25,  3.  ],
            dtype=float32)
    • time
      (time)
      datetime64[ns]
      2018-01-01T23:00:00 ... 2022-12-...
      long_name :
      time
      array(['2018-01-01T23:00:00.000000000', '2018-02-01T23:00:00.000000000',
             '2018-03-01T23:00:00.000000000', '2018-04-01T23:00:00.000000000',
             '2018-05-01T23:00:00.000000000', '2018-06-01T23:00:00.000000000',
             '2018-07-01T23:00:00.000000000', '2018-08-01T23:00:00.000000000',
             '2018-09-01T23:00:00.000000000', '2018-10-01T23:00:00.000000000',
             '2018-11-01T23:00:00.000000000', '2018-12-01T23:00:00.000000000',
             '2019-01-01T23:00:00.000000000', '2019-02-01T23:00:00.000000000',
             '2019-03-01T23:00:00.000000000', '2019-04-01T23:00:00.000000000',
             '2019-05-01T23:00:00.000000000', '2019-06-01T23:00:00.000000000',
             '2019-07-01T23:00:00.000000000', '2019-08-01T23:00:00.000000000',
             '2019-09-01T23:00:00.000000000', '2019-10-01T23:00:00.000000000',
             '2019-11-01T23:00:00.000000000', '2019-12-01T23:00:00.000000000',
             '2020-01-01T23:00:00.000000000', '2020-02-01T23:00:00.000000000',
             '2020-03-01T23:00:00.000000000', '2020-04-01T23:00:00.000000000',
             '2020-05-01T23:00:00.000000000', '2020-06-01T23:00:00.000000000',
             '2020-07-01T23:00:00.000000000', '2020-08-01T23:00:00.000000000',
             '2020-09-01T23:00:00.000000000', '2020-10-01T23:00:00.000000000',
             '2020-11-01T23:00:00.000000000', '2020-12-01T23:00:00.000000000',
             '2021-01-01T23:00:00.000000000', '2021-02-01T23:00:00.000000000',
             '2021-03-01T23:00:00.000000000', '2021-04-01T23:00:00.000000000',
             '2021-05-01T23:00:00.000000000', '2021-06-01T23:00:00.000000000',
             '2021-07-01T23:00:00.000000000', '2021-08-01T23:00:00.000000000',
             '2021-09-01T23:00:00.000000000', '2021-10-01T23:00:00.000000000',
             '2021-11-01T23:00:00.000000000', '2021-12-01T23:00:00.000000000',
             '2022-01-01T23:00:00.000000000', '2022-02-01T23:00:00.000000000',
             '2022-03-01T23:00:00.000000000', '2022-04-01T23:00:00.000000000',
             '2022-05-01T23:00:00.000000000', '2022-06-01T23:00:00.000000000',
             '2022-07-01T23:00:00.000000000', '2022-08-01T23:00:00.000000000',
             '2022-09-01T23:00:00.000000000', '2022-10-01T23:00:00.000000000',
             '2022-11-01T23:00:00.000000000', '2022-12-01T23:00:00.000000000'],
            dtype='datetime64[ns]')
    • t2m
      (time, latitude, longitude)
      float32
      ...
      units :
      K
      long_name :
      2 metre temperature
      [179340 values with dtype=float32]
    • sst
      (time, latitude, longitude)
      float32
      ...
      units :
      K
      long_name :
      Sea surface temperature
      [179340 values with dtype=float32]
    • tp
      (time, latitude, longitude)
      float32
      ...
      units :
      m
      long_name :
      Total precipitation
      [179340 values with dtype=float32]
  • Conventions :
    CF-1.6
    history :
    2023-03-14 04:44:29 GMT by grib_to_netcdf-2.25.1: /opt/ecmwf/mars-client/bin/grib_to_netcdf.bin -S param -o /cache/data4/adaptor.mars.internal-1678769068.6900754-32437-17-b23b0f81-0686-4bdf-b4ef-3e356f6f67de.nc /cache/tmp/b23b0f81-0686-4bdf-b4ef-3e356f6f67de-adaptor.mars.internal-1678769038.9013975-32437-26-tmp.grib
In [45]:
plt.figure(figsize=(12, 8))

extent = [33, 48, 3, 15]


lon2d, lat2d = np.meshgrid(era5_ecmwf_mon_ds.longitude,era5_ecmwf_mon_ds.latitude)

ax11 = plt.axes(projection=ccrs.PlateCarree())

    
c1= ax11.pcolormesh(lon2d, lat2d,era5_ecmwf_mon_ds.t2m[0,:,:], 
                    transform=ccrs.PlateCarree(), cmap='viridis_r')

ax11.coastlines()

ax11.add_feature(cfeature.OCEAN, zorder=100, edgecolor='k')
ax11.add_feature(cfeature.BORDERS, edgecolor='black')
ax11.add_feature(cfeature.LAKES, edgecolor='black')


ax11.set_extent(extent)
ax11.gridlines()

#grid_lines = ax11.gridlines(draw_labels=True)

#grid_lines.xformater= LONGITUDE_FORMATTER
#grid_lines.yformater= LATITUDE_FORMATTER

ax11.set_xticks(np.arange(33,48,2), crs=ccrs.PlateCarree())
ax11.set_yticks(np.arange(3,15,2), crs=ccrs.PlateCarree())

lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()

ax11.xaxis.set_major_formatter(lon_formatter)
ax11.yaxis.set_major_formatter(lat_formatter)

ax11.set_xlabel('Longitude [degrees east]', fontsize=12)
ax11.set_ylabel('Latitude [degrees north]', fontsize=12)

cb = plt.colorbar(c1,orientation="vertical",extendrect='True', shrink=0.90)
cb.set_label("temperature [K]", fontsize=12)

plt.title("2 metre temperature", fontsize=16, pad=15.0)
plt.subplots_adjust(top=0.92)

#plt.savefig('2metretemperature.jpeg', transparent=True,  bbox_inches='tight', dpi=800)

plt.show()
#contourf
#pcolor
#contour
#pcolormesh

Sub plot_1¶

In [46]:
fg = era5_ecmwf_mon_ds.t2m.isel(time=slice(0, 12)).plot(col='time', col_wrap=4,
        transform=ccrs.PlateCarree(), subplot_kws=dict(projection=ccrs.PlateCarree()), 
                        figsize=(20, 18),  cmap='cool',  cbar_kwargs={'shrink': 0.7})
                 
for ax in fg.axes.flat:
    ax.set_extent([33, 48, 3, 15], ccrs.PlateCarree())
    ax.coastlines()
    ax.add_feature(cfeature.OCEAN, zorder=100, edgecolor='k', facecolor='blue')
    ax.add_feature(cfeature.BORDERS, edgecolor='black')
    ax.add_feature(cfeature.LAKES, edgecolor='black') 

Sub plot 2¶

In [47]:
plt.figure(figsize=(12, 8))

plt.rcParams['axes.titlesize'] =17
plt.rcParams['axes.titlepad'] =10

plt.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0.3, hspace=0.4)

plt.suptitle("Sub Plot ",  fontsize=24)

lon2d, lat2d = np.meshgrid(era5_ecmwf_mon_ds.longitude,era5_ecmwf_mon_ds.latitude)  

extent = [33, 48, 3, 15]
central_lon, central_lat = 37,7

ax1 = plt.subplot(1,4,1, projection=ccrs.PlateCarree(central_longitude=central_lon, globe=None))
#Bare_soil.plot.pcolormesh(ax=ax1, transform=ccrs.PlateCarree(), cmap="tab20c", vmin=0, vmax=100)
ct1 = ax1.pcolormesh(lon2d, lat2d, era5_ecmwf_mon_ds.t2m[0,:,:], 
                     transform=ccrs.PlateCarree(), cmap='tab20c')
ax1.add_feature(cfeature.OCEAN, zorder=100, edgecolor='k', facecolor='blue')
ax1.add_feature(cfeature.BORDERS, edgecolor='black')
ax1.add_feature(cfeature.LAKES, edgecolor='black')
ax1.coastlines(resolution='10m')
ax1.set_extent(extent)
ax1.gridlines()
ax1.set_xticks(np.arange(33,49,3), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(3,16,3), crs=ccrs.PlateCarree())

lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
#cb = plt.colorbar(ct1 ,orientation="vertical",extendrect='True', shrink = 0.1)
#cb.set_label("K", fontsize=12)
ax1.set(title="Month 1")
        
ax2 = plt.subplot(1,4,2, projection=ccrs.PlateCarree(central_longitude=central_lon, globe=None))
#NET_Temperate.plot.pcolormesh(ax=ax2,transform=ccrs.PlateCarree(), cmap="tab20c", vmin=0, vmax=100)
ct2 = ax2.pcolormesh(lon2d, lat2d, era5_ecmwf_mon_ds.t2m[1,:,:],
                     transform=ccrs.PlateCarree(), cmap='tab20c')
ax2.add_feature(cfeature.OCEAN, zorder=100, edgecolor='k', facecolor='blue')
ax2.add_feature(cfeature.BORDERS, edgecolor='black')
ax2.add_feature(cfeature.LAKES, edgecolor='black')
ax2.coastlines(resolution='10m')
ax2.set_extent(extent)
ax2.gridlines()
ax2.set_xticks(np.arange(33,49,3), crs=ccrs.PlateCarree())
ax2.set_yticks(np.arange(3,16,3), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax2.xaxis.set_major_formatter(lon_formatter)
ax2.yaxis.set_major_formatter(lat_formatter)
cb = plt.colorbar(ct2 ,orientation="vertical",extendrect='True', shrink = 0.1)
#cb.set_label("K", fontsize=12)
ax2.set(title="Month 2")

ax3 = plt.subplot(1,4,3, projection=ccrs.PlateCarree(central_longitude=central_lon, globe=None))
#NET_Boreal.plot.pcolormesh(ax=ax3,transform=ccrs.PlateCarree(), cmap="tab20c", vmin=0, vmax=100)
ct3 = ax3.pcolormesh(lon2d, lat2d, era5_ecmwf_mon_ds.t2m[2,:,:], 
                     transform=ccrs.PlateCarree(), cmap='tab20c')
ax3.add_feature(cfeature.OCEAN, zorder=100, edgecolor='k', facecolor='blue')
ax3.add_feature(cfeature.BORDERS, edgecolor='black')
ax3.add_feature(cfeature.LAKES, edgecolor='black')
ax3.coastlines(resolution='10m')
ax3.set_extent(extent)
ax3.gridlines()
ax3.set_xticks(np.arange(33,49,3), crs=ccrs.PlateCarree())
ax3.set_yticks(np.arange(3,16,3), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax3.xaxis.set_major_formatter(lon_formatter)
ax3.yaxis.set_major_formatter(lat_formatter)
cb = plt.colorbar(ct3 ,orientation="vertical",extendrect='True', shrink = 0.1)
# cb.set_label("K", fontsize=12)
ax3.set(title="Month 3")


ax4 = plt.subplot(1,4,4, projection=ccrs.PlateCarree(central_longitude=central_lon, globe=None))
#NDT_Boreal.plot.pcolormesh(ax=ax4,transform=ccrs.PlateCarree(), cmap="tab20c", vmin=0, vmax=100)
ct4 = ax4.pcolormesh(lon2d, lat2d, era5_ecmwf_mon_ds.t2m[3,:,:], 
                     transform=ccrs.PlateCarree(), cmap='tab20c')
ax4.add_feature(cfeature.OCEAN, zorder=100, edgecolor='k', facecolor='blue')
ax4.add_feature(cfeature.BORDERS, edgecolor='black')
ax4.add_feature(cfeature.LAKES, edgecolor='black')
ax4.coastlines(resolution='10m')
ax4.set_extent(extent)
ax4.gridlines()
ax4.set_xticks(np.arange(33,49,3), crs=ccrs.PlateCarree())
ax4.set_yticks(np.arange(3,16,3), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax4.xaxis.set_major_formatter(lon_formatter)
ax4.yaxis.set_major_formatter(lat_formatter)
cb = plt.colorbar(ct4 ,orientation="vertical",extendrect='True', shrink = 0.1)
cb.set_label("K", fontsize=12 )
ax4.set(title="Month 4")

plt.subplots_adjust(top=1.45)
#plt.savefig('sub.jpeg', transparent=True,  bbox_inches='tight', dpi=800)

plt.show()

Color option¶

['Perceptually Uniform Sequential'] = ['viridis', 'plasma', 'inferno', 'magma', 'cividis']

['Sequential'] = ['Greys', 'Purples', 'Blues', 'Greens', 'Oranges', 'Reds','YlOrBr',

'YlOrRd', 'OrRd', 'PuRd', 'RdPu', 'BuPu','GnBu', 'PuBu', 'YlGnBu', 'PuBuGn', 'BuGn', 'YlGn']¶

['Perceptually Uniform Sequential'] = ['viridis', 'plasma', 'inferno', 'magma', 'cividis']

['Sequential'] = ['Greys', 'Purples', 'Blues', 'Greens', 'Oranges', 'Reds','YlOrBr', 'YlOrRd','OrRd', 'PuRd', 'RdPu', 'BuPu', 'GnBu', 'PuBu', 'YlGnBu', 'PuBuGn', 'BuGn', 'YlGn']

['Sequential (2)'] = ['binary', 'gist_yarg', 'gist_gray', 'gray', 'bone', 'pink', 'spring', 'summer', 'autumn', 'winter', 'cool', 'Wistia','hot', 'afmhot', 'gist_heat', 'copper']

['Diverging'] = ['PiYG', 'PRGn', 'BrBG', 'PuOr', 'RdGy', 'RdBu','RdYlBu', 'RdYlGn', 'Spectral', 'coolwarm', 'bwr', 'seismic']

['Cyclic'] = ['twilight', 'twilight_shifted', 'hsv']['Qualitative'] = ['Pastel1', 'Pastel2', 'Paired', 'Accent','Dark2', 'Set1', 'Set2', 'Set3', 'tab20c', 'tab20', 'tab20b', 'tab20c']

['Miscellaneous'] = ['flag', 'prism', 'ocean', 'gist_earth', 'terrain', 'gist_stern','gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'brg','gist_rainbow', 'rainbow', 'jet', 'nipy_spectral', 'gist_ncar']

In [ ]: